Making statistics work for (evolutionary language) scientists

TTF DataScience, WP Statistics

Erik J. Ringen

Introduction

This session is about how statistics can be in service of science.

  • Key ideas:
    • Scientist learn the language of statistical models in order to do their work
    • But a lot gets lost in translation when we try to translate back to substantive scientific questions from the output of statistical models
  • Goals:
    • Part 1: Teach you a general method to translate statistical models back into substantive scientific questions using marginaleffects
    • Part 2: Managing the workflow of statistical computing using targets

What this session is not about

  • Promoting any particular statistical model
  • Telling you what you ‘should’ be studying (statistician’s fallacy)
  • Bayesian statistics

Software

In this session I will use:

  • R (latest version)
  • The following R packages:
    • marginaleffects
    • tidyverse
    • targets
    • mgcv (tangentially)
    • knitr (for rendering outputs)

Quarto markdown and code here (optional):

Part 1: Marginal effects

About me

  • I earned my PhD in Evolutionary Anthropology during the 433rd lunation of the third millennium, when Mars was at an orbital longitude of ~311 degrees in May 2022

  • To get from Zurich to here, I travelled 807 light-microseconds 242 kilometers

  • The units we use to communicate matter!

The way we report statistics is often confusing (and misleading)

The way we report statistics is often confusing (and misleading)

In contrast, descriptive statistics and data viz are easy to understand

‘Monolingual participants have a mean cognitive flexibility score of 51.27 (SD = 18.13), while multilingual participants have a mean score of 61.77 (SD = 18.78).’

  • Simple comparisons with scientifically-meaningful units.
  • What if we could do the same for statistical models?

Predictions, not parameters

  • The central idea: rather than directly interpret parameters, use the model to make predictions
    • contrasts between different predictions consitute marginal effects
  • The marginaleffects package by Vincent Arel-Bundock makes these calculations easy for a wide range of models in both R and Python.

In simple models, a parameter might be all we need

  • Example: do multilingual children have greater cognitive flexibility than monolingual children?
    • Data: 200 children from some population
    • Estimand: what is the average difference in cognitive flexibility (measured by a psychometric instrument ranging from 0 to 100) between multilingual and monolingual children?
    • Estimate: ‘Multilingual participants have an average cognitive flexibility score of 10.508 [5.907, 15.108] points higher than monolingual participants.’


Call:
lm(formula = cognitive_flexibility ~ multilingual)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.387 -11.913   1.581  13.388  36.985 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)    51.267      1.665  30.795  < 2e-16 ***
multilingual   10.508      2.336   4.499 1.05e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.46 on 248 degrees of freedom
Multiple R-squared:  0.07545,   Adjusted R-squared:  0.07172 
F-statistic: 20.24 on 1 and 248 DF,  p-value: 1.051e-05

A regression is a machine that makes predictions

\[ \text{Cognitive flexibility} = \alpha + \beta \times \text{Multilingual} + \epsilon \] \[ \epsilon \sim \text{Normal}(0, \sigma^2) \]

  • The parts of the machine are the parameters (\(\alpha, \beta, \sigma\)):
    • To make a prediction for a monolingual child, just take the intercept (\(\alpha\)).
    • To make a prediction for a multilingual child, add the coefficient for multilingual (\(\beta\)) to the intercept.
  • Linear regression is such a simple machine that it offers a 1:1 mapping between the parameters and the estimand (average difference between multilingual and monolingual children).

Interactions

Interactions imply heterogeneity

  • For a given child, let \(Y_i(\text{multilingual})\) be the outcome if the child is multilingual and \(Y_i(\text{monolingual})\) be the outcome if the child is monolingual.

  • The estimand is the average difference in potential outcomes: \(\frac{1}{N} \sum_{i=1}^N [Y_i(\text{multilingual}) - Y_i(\text{monolingual})]\). Known as the average treatment effect (ATE), or average marginal effect (AME).

  • In our simple linear regression, all children are identical, conditional on their mono/multilingual status, which allows us to reduce the estimand to a single parameter.

Interactions imply heterogeneity


Call:
lm(formula = cognitive_flexibility ~ multilingual * age_c)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.937  -9.753   1.457  12.291  38.200 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)          51.123      1.580  32.353  < 2e-16 ***
multilingual         10.969      2.217   4.948 1.39e-06 ***
age_c                 2.316      1.079   2.145   0.0329 *  
multilingual:age_c    2.949      1.506   1.959   0.0512 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.51 on 246 degrees of freedom
Multiple R-squared:  0.1752,    Adjusted R-squared:  0.1652 
F-statistic: 17.42 on 3 and 246 DF,  p-value: 2.738e-10

  • Which coefficient corresponds to the estimand (average marginal effect)?
    • None! The AME depends on the distribution of the moderator (age):
      • The magnitude of the multilingual effect depends on age–now it matters which children we are averaging over.

Meaning of interaction term is highly contextual

On the previous slide, the interaction term was mean-centered. What if we use the raw age variable instead?


Call:
lm(formula = cognitive_flexibility ~ multilingual * age)

Residuals:
    Min      1Q  Median      3Q     Max 
-45.937  -9.753   1.457  12.291  38.200 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)   
(Intercept)        29.321     10.350   2.833   0.0050 **
multilingual      -16.799     14.351  -1.171   0.2429   
age                 2.316      1.079   2.145   0.0329 * 
multilingual:age    2.949      1.506   1.959   0.0512 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.51 on 246 degrees of freedom
Multiple R-squared:  0.1752,    Adjusted R-squared:  0.1652 
F-statistic: 17.42 on 3 and 246 DF,  p-value: 2.738e-10

Getting started with marginaleffects: predictions

library(marginaleffects)

uncentered_model <- lm(cognitive_flexibility ~ multilingual * age, data = d)

preds <- predictions(uncentered_model)

preds |> 
    as.data.frame() |> 
    dplyr::select(rowid, multilingual, age, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)
rowid multilingual age estimate conf.low conf.high
1 monolingual 8.29 48.51 44.52 52.50
2 monolingual 8.11 48.10 43.86 52.34
3 monolingual 9.97 52.40 49.14 55.66
4 multilingual 8.34 56.42 52.73 60.12
5 monolingual 9.66 51.68 48.56 54.80
6 multilingual 10.93 70.05 65.61 74.49
7 multilingual 7.84 53.80 49.45 58.16
8 multilingual 9.02 60.02 56.90 63.15
9 monolingual 9.36 50.99 47.89 54.10
10 monolingual 11.34 55.58 50.57 60.59
11 multilingual 11.63 73.75 68.17 79.33
12 multilingual 11.41 72.60 67.39 77.81
13 multilingual 10.37 67.13 63.43 70.82
14 monolingual 11.75 56.53 50.81 62.25
15 multilingual 9.58 62.97 59.89 66.06
16 monolingual 9.88 52.21 49.00 55.42
17 multilingual 8.68 58.23 54.89 61.58
18 monolingual 8.74 49.55 46.08 53.02
19 monolingual 7.10 45.76 39.86 51.67
20 monolingual 9.51 51.35 48.26 54.45
21 monolingual 11.36 55.62 50.58 60.65
22 multilingual 7.03 49.54 43.88 55.21
23 monolingual 7.36 46.37 40.92 51.81
24 monolingual 7.82 47.43 42.76 52.11
25 monolingual 10.85 54.45 50.20 58.70
26 monolingual 10.68 54.04 50.04 58.04
27 multilingual 11.86 74.96 68.98 80.95
28 multilingual 9.33 61.66 58.61 64.70
29 monolingual 7.37 46.39 40.97 51.82
30 multilingual 10.24 66.46 62.91 70.01
31 monolingual 10.79 54.31 50.15 58.48
32 multilingual 7.69 52.99 48.40 57.58
33 monolingual 8.98 50.12 46.86 53.39
34 multilingual 8.12 55.30 51.34 59.26
35 multilingual 7.29 50.90 45.68 56.13
36 monolingual 8.98 50.12 46.85 53.38
37 monolingual 7.32 46.28 40.78 51.79
38 monolingual 8.13 48.15 43.94 52.35
39 monolingual 7.27 46.16 40.57 51.76
40 multilingual 10.35 67.02 63.35 70.70
41 monolingual 8.49 48.98 45.24 52.71
42 multilingual 7.50 52.03 47.15 56.91
43 multilingual 7.36 51.27 46.16 56.38
44 monolingual 11.40 55.73 50.61 60.84
45 monolingual 10.77 54.26 50.13 58.40
46 monolingual 11.08 54.99 50.39 59.58
47 monolingual 11.91 56.90 50.90 62.91
48 multilingual 7.52 52.11 47.25 56.96
49 monolingual 7.50 46.68 41.47 51.89
50 monolingual 10.99 54.78 50.32 59.24
51 multilingual 10.92 70.03 65.60 74.47
52 monolingual 7.05 45.64 39.64 51.64
53 monolingual 10.90 54.55 50.24 58.86
54 monolingual 10.65 53.98 50.01 57.94
55 monolingual 10.15 52.83 49.42 56.23
56 multilingual 9.40 62.04 58.99 65.09
57 multilingual 7.78 53.50 49.06 57.94
58 monolingual 7.04 45.63 39.62 51.64
59 multilingual 9.26 61.29 58.24 64.34
60 monolingual 9.46 51.23 48.14 54.33
61 monolingual 8.95 50.04 46.75 53.33
62 multilingual 9.32 61.61 58.57 64.66
63 multilingual 10.57 68.16 64.22 72.09
64 monolingual 7.28 46.17 40.58 51.76
65 monolingual 8.77 49.64 46.21 53.07
66 multilingual 11.01 70.51 65.94 75.09
67 monolingual 11.18 55.21 50.46 59.95
68 monolingual 8.19 48.28 44.16 52.41
69 monolingual 8.77 49.63 46.19 53.07
70 monolingual 11.28 55.45 50.53 60.37
71 multilingual 11.27 71.85 66.88 76.83
72 monolingual 8.48 48.96 45.21 52.70
73 monolingual 7.74 47.23 42.42 52.05
74 monolingual 10.52 53.68 49.88 57.48
75 monolingual 7.52 46.73 41.56 51.90
76 multilingual 7.17 50.27 44.84 55.70
77 multilingual 12.00 75.69 69.46 81.92
78 monolingual 7.17 45.94 40.16 51.71
79 multilingual 8.69 58.29 54.95 61.62
80 multilingual 11.58 73.47 67.98 78.96
81 multilingual 10.09 65.63 62.23 69.02
82 multilingual 8.43 56.92 53.33 60.50
83 monolingual 10.69 54.07 50.06 58.09
84 multilingual 11.17 71.34 66.52 76.16
85 multilingual 8.57 57.65 54.21 61.10
86 multilingual 9.46 62.35 59.29 65.40
87 monolingual 10.49 53.61 49.85 57.37
88 monolingual 10.21 52.96 49.50 56.42
89 multilingual 10.22 66.33 62.80 69.86
90 monolingual 11.89 56.85 50.89 62.82
91 multilingual 9.07 60.30 57.20 63.40
92 multilingual 7.60 52.52 47.79 57.25
93 monolingual 9.63 51.62 48.51 54.73
94 multilingual 8.13 55.30 51.35 59.26
95 multilingual 9.43 62.18 59.13 65.23
96 monolingual 8.85 49.82 46.45 53.18
97 monolingual 11.92 56.92 50.90 62.94
98 multilingual 8.94 59.60 56.44 62.76
99 monolingual 8.15 48.19 44.00 52.37
100 monolingual 10.12 52.75 49.37 56.13
101 monolingual 7.68 47.11 42.22 52.01
102 multilingual 11.84 74.85 68.90 80.80
103 multilingual 9.58 62.94 59.86 66.02
104 monolingual 7.82 47.42 42.74 52.10
105 multilingual 10.11 65.75 62.33 69.17
106 multilingual 11.93 75.33 69.22 81.45
107 monolingual 10.34 53.27 49.68 56.87
108 monolingual 9.09 50.38 47.18 53.58
109 monolingual 8.62 49.28 45.69 52.86
110 multilingual 11.18 71.37 66.54 76.20
111 multilingual 7.72 53.16 48.63 57.70
112 multilingual 7.96 54.45 50.28 58.63
113 multilingual 11.48 72.99 67.65 78.32
114 monolingual 8.54 49.10 45.42 52.77
115 monolingual 8.82 49.74 46.34 53.13
116 multilingual 10.92 70.02 65.58 74.45
117 multilingual 7.97 54.47 50.30 58.64
118 multilingual 7.09 49.85 44.28 55.41
119 monolingual 9.03 50.24 47.01 53.47
120 multilingual 9.42 62.10 59.05 65.15
121 monolingual 9.11 50.42 47.23 53.61
122 multilingual 8.71 58.40 55.09 61.72
123 monolingual 11.33 55.56 50.57 60.56
124 monolingual 9.28 50.80 47.68 53.92
125 monolingual 9.67 51.71 48.59 54.83
126 multilingual 11.82 74.75 68.84 80.67
127 multilingual 10.87 69.77 65.41 74.13
128 monolingual 8.04 47.95 43.62 52.28
129 monolingual 8.54 49.11 45.44 52.78
130 monolingual 11.86 56.78 50.87 62.69
131 multilingual 9.92 64.78 61.51 68.04
132 monolingual 10.80 54.34 50.16 58.52
133 multilingual 8.86 59.19 55.98 62.40
134 multilingual 10.85 69.63 65.31 73.95
135 multilingual 9.69 63.53 60.41 66.66
136 multilingual 11.57 73.44 67.96 78.92
137 multilingual 7.93 54.26 50.03 58.49
138 multilingual 8.41 56.81 53.20 60.42
139 multilingual 7.47 51.88 46.96 56.80
140 monolingual 8.05 47.97 43.65 52.29
141 multilingual 11.89 75.10 69.07 81.13
142 multilingual 8.48 57.18 53.64 60.71
143 multilingual 10.63 68.49 64.47 72.51
144 monolingual 10.93 54.63 50.27 58.99
145 multilingual 7.53 52.15 47.32 56.99
146 multilingual 8.20 55.69 51.82 59.55
147 multilingual 8.35 56.50 52.82 60.18
148 monolingual 7.51 46.70 41.51 51.89
149 multilingual 7.59 52.48 47.74 57.22
150 multilingual 11.96 75.47 69.32 81.63
151 monolingual 11.93 56.95 50.91 62.99
152 multilingual 7.69 52.99 48.40 57.58
153 monolingual 11.53 56.01 50.69 61.34
154 multilingual 9.88 64.55 61.32 67.78
155 monolingual 8.98 50.11 46.84 53.38
156 monolingual 9.25 50.74 47.61 53.87
157 monolingual 10.53 53.71 49.90 57.53
158 monolingual 7.41 46.49 41.13 51.84
159 monolingual 8.70 49.46 45.95 52.97
160 monolingual 10.40 53.41 49.75 57.08
161 monolingual 8.58 49.20 45.58 52.83
162 multilingual 11.16 71.27 66.47 76.07
163 multilingual 8.08 55.04 51.02 59.07
164 multilingual 9.49 62.49 59.43 65.55
165 monolingual 8.38 48.73 44.86 52.59
166 multilingual 7.96 54.43 50.25 58.62
167 monolingual 11.75 56.54 50.81 62.26
168 monolingual 8.61 49.26 45.66 52.85
169 multilingual 9.39 61.97 58.93 65.02
170 monolingual 7.14 45.86 40.02 51.69
171 multilingual 9.74 63.79 60.65 66.94
172 multilingual 10.22 66.34 62.81 69.87
173 multilingual 9.98 65.08 61.77 68.38
174 multilingual 8.61 57.85 54.44 61.26
175 monolingual 11.46 55.85 50.64 61.05
176 monolingual 10.13 52.78 49.39 56.17
177 monolingual 8.51 49.04 45.33 52.74
178 multilingual 8.94 59.60 56.44 62.76
179 multilingual 7.80 53.60 49.19 58.02
180 multilingual 11.31 72.09 67.04 77.14
181 multilingual 11.77 74.47 68.65 80.29
182 monolingual 9.82 52.06 48.88 55.24
183 monolingual 8.65 49.35 45.79 52.90
184 monolingual 11.98 57.07 50.93 63.21
185 monolingual 8.17 48.25 44.11 52.39
186 multilingual 10.06 65.51 62.13 68.88
187 multilingual 7.54 52.23 47.41 57.04
188 monolingual 9.44 51.17 48.07 54.27
189 monolingual 7.50 46.68 41.48 51.89
190 multilingual 7.81 53.62 49.21 58.03
191 multilingual 8.41 56.83 53.22 60.44
192 monolingual 9.92 52.29 49.06 55.52
193 monolingual 10.66 54.00 50.03 57.98
194 multilingual 7.83 53.74 49.36 58.11
195 multilingual 11.33 72.19 67.11 77.27
196 monolingual 10.54 53.74 49.91 57.56
197 monolingual 10.80 54.34 50.16 58.51
198 multilingual 7.74 53.25 48.74 57.76
199 monolingual 8.79 49.68 46.26 53.10
200 multilingual 10.37 67.10 63.42 70.79
201 multilingual 9.62 63.17 60.08 66.26
202 multilingual 8.75 58.59 55.30 61.88
203 multilingual 8.20 55.71 51.85 59.57
204 monolingual 7.29 46.21 40.64 51.77
205 multilingual 8.18 55.61 51.72 59.49
206 multilingual 11.45 72.81 67.53 78.09
207 multilingual 11.06 70.75 66.11 75.40
208 monolingual 10.74 54.19 50.10 58.27
209 multilingual 7.77 53.46 49.00 57.91
210 multilingual 7.62 52.66 47.98 57.35
211 multilingual 11.87 75.04 69.03 81.05
212 monolingual 9.18 50.58 47.42 53.74
213 monolingual 9.32 50.90 47.79 54.02
214 multilingual 7.83 53.73 49.35 58.11
215 multilingual 9.92 64.78 61.51 68.04
216 multilingual 8.35 56.51 52.83 60.18
217 multilingual 8.15 55.44 51.51 59.36
218 monolingual 10.46 53.53 49.81 57.26
219 multilingual 8.41 56.82 53.22 60.43
220 monolingual 11.05 54.91 50.37 59.46
221 multilingual 7.47 51.85 46.92 56.78
222 monolingual 11.11 55.05 50.41 59.69
223 monolingual 9.14 50.48 47.30 53.66
224 monolingual 10.78 54.28 50.14 58.43
225 multilingual 10.31 66.82 63.19 70.44
226 multilingual 9.22 61.08 58.02 64.14
227 monolingual 10.14 52.79 49.40 56.19
228 monolingual 7.00 45.54 39.46 51.62
229 multilingual 8.09 55.10 51.09 59.11
230 monolingual 10.52 53.69 49.89 57.50
231 monolingual 8.08 48.02 43.74 52.31
232 multilingual 11.07 70.81 66.14 75.47
233 multilingual 8.54 57.48 54.00 60.96
234 monolingual 10.44 53.49 49.79 57.20
235 multilingual 11.66 73.93 68.29 79.58
236 multilingual 7.58 52.43 47.67 57.18
237 multilingual 7.64 52.74 48.08 57.40
238 multilingual 10.39 67.23 63.52 70.95
239 multilingual 9.14 60.67 57.60 63.75
240 monolingual 11.17 55.19 50.46 59.93
241 monolingual 11.86 56.78 50.87 62.69
242 multilingual 7.35 51.23 46.11 56.36
243 monolingual 9.30 50.85 47.74 53.97
244 multilingual 10.51 67.85 63.99 71.71
245 monolingual 7.43 46.54 41.22 51.85
246 monolingual 11.96 57.03 50.92 63.13
247 monolingual 8.27 48.46 44.44 52.48
248 multilingual 7.25 50.68 45.39 55.98
249 multilingual 10.43 67.45 63.68 71.21
250 multilingual 10.93 70.10 65.64 74.55

Getting started with marginaleffects: predictions

library(marginaleffects)

centered_model <- lm(cognitive_flexibility ~ multilingual * age_c, data = d)

preds <- predictions(centered_model)

preds |> 
    as.data.frame() |> 
    dplyr::select(rowid, multilingual, age_c, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)
rowid multilingual age_c estimate conf.low conf.high
1 monolingual -1.13 48.51 44.52 52.50
2 monolingual -1.31 48.10 43.86 52.34
3 monolingual 0.55 52.40 49.14 55.66
4 multilingual -1.08 56.42 52.73 60.12
5 monolingual 0.24 51.68 48.56 54.80
6 multilingual 1.51 70.05 65.61 74.49
7 multilingual -1.57 53.80 49.45 58.16
8 multilingual -0.39 60.02 56.90 63.15
9 monolingual -0.06 50.99 47.89 54.10
10 monolingual 1.93 55.58 50.57 60.59
11 multilingual 2.21 73.75 68.17 79.33
12 multilingual 2.00 72.60 67.39 77.81
13 multilingual 0.96 67.13 63.43 70.82
14 monolingual 2.34 56.53 50.81 62.25
15 multilingual 0.17 62.97 59.89 66.06
16 monolingual 0.47 52.21 49.00 55.42
17 multilingual -0.73 58.23 54.89 61.58
18 monolingual -0.68 49.55 46.08 53.02
19 monolingual -2.31 45.76 39.86 51.67
20 monolingual 0.10 51.35 48.26 54.45
21 monolingual 1.94 55.62 50.58 60.65
22 multilingual -2.38 49.54 43.88 55.21
23 monolingual -2.05 46.37 40.92 51.81
24 monolingual -1.59 47.43 42.76 52.11
25 monolingual 1.44 54.45 50.20 58.70
26 monolingual 1.26 54.04 50.04 58.04
27 multilingual 2.44 74.96 68.98 80.95
28 multilingual -0.08 61.66 58.61 64.70
29 monolingual -2.04 46.39 40.97 51.82
30 multilingual 0.83 66.46 62.91 70.01
31 monolingual 1.38 54.31 50.15 58.48
32 multilingual -1.73 52.99 48.40 57.58
33 monolingual -0.43 50.12 46.86 53.39
34 multilingual -1.29 55.30 51.34 59.26
35 multilingual -2.12 50.90 45.68 56.13
36 monolingual -0.44 50.12 46.85 53.38
37 monolingual -2.09 46.28 40.78 51.79
38 monolingual -1.29 48.15 43.94 52.35
39 monolingual -2.14 46.16 40.57 51.76
40 multilingual 0.94 67.02 63.35 70.70
41 monolingual -0.93 48.98 45.24 52.71
42 multilingual -1.91 52.03 47.15 56.91
43 multilingual -2.06 51.27 46.16 56.38
44 monolingual 1.99 55.73 50.61 60.84
45 monolingual 1.36 54.26 50.13 58.40
46 monolingual 1.67 54.99 50.39 59.58
47 monolingual 2.50 56.90 50.90 62.91
48 multilingual -1.90 52.11 47.25 56.96
49 monolingual -1.92 46.68 41.47 51.89
50 monolingual 1.58 54.78 50.32 59.24
51 multilingual 1.51 70.03 65.60 74.47
52 monolingual -2.37 45.64 39.64 51.64
53 monolingual 1.48 54.55 50.24 58.86
54 monolingual 1.23 53.98 50.01 57.94
55 monolingual 0.74 52.83 49.42 56.23
56 multilingual -0.01 62.04 58.99 65.09
57 multilingual -1.63 53.50 49.06 57.94
58 monolingual -2.37 45.63 39.62 51.64
59 multilingual -0.15 61.29 58.24 64.34
60 monolingual 0.05 51.23 48.14 54.33
61 monolingual -0.47 50.04 46.75 53.33
62 multilingual -0.09 61.61 58.57 64.66
63 multilingual 1.15 68.16 64.22 72.09
64 monolingual -2.14 46.17 40.58 51.76
65 monolingual -0.64 49.64 46.21 53.07
66 multilingual 1.60 70.51 65.94 75.09
67 monolingual 1.76 55.21 50.46 59.95
68 monolingual -1.23 48.28 44.16 52.41
69 monolingual -0.64 49.63 46.19 53.07
70 monolingual 1.87 55.45 50.53 60.37
71 multilingual 1.85 71.85 66.88 76.83
72 monolingual -0.94 48.96 45.21 52.70
73 monolingual -1.68 47.23 42.42 52.05
74 monolingual 1.11 53.68 49.88 57.48
75 monolingual -1.90 46.73 41.56 51.90
76 multilingual -2.25 50.27 44.84 55.70
77 multilingual 2.58 75.69 69.46 81.92
78 monolingual -2.24 45.94 40.16 51.71
79 multilingual -0.72 58.29 54.95 61.62
80 multilingual 2.16 73.47 67.98 78.96
81 multilingual 0.67 65.63 62.23 69.02
82 multilingual -0.98 56.92 53.33 60.50
83 monolingual 1.27 54.07 50.06 58.09
84 multilingual 1.76 71.34 66.52 76.16
85 multilingual -0.84 57.65 54.21 61.10
86 multilingual 0.05 62.35 59.29 65.40
87 monolingual 1.07 53.61 49.85 57.37
88 monolingual 0.79 52.96 49.50 56.42
89 multilingual 0.81 66.33 62.80 69.86
90 monolingual 2.47 56.85 50.89 62.82
91 multilingual -0.34 60.30 57.20 63.40
92 multilingual -1.82 52.52 47.79 57.25
93 monolingual 0.22 51.62 48.51 54.73
94 multilingual -1.29 55.30 51.35 59.26
95 multilingual 0.02 62.18 59.13 65.23
96 monolingual -0.56 49.82 46.45 53.18
97 monolingual 2.50 56.92 50.90 62.94
98 multilingual -0.47 59.60 56.44 62.76
99 monolingual -1.27 48.19 44.00 52.37
100 monolingual 0.70 52.75 49.37 56.13
101 monolingual -1.73 47.11 42.22 52.01
102 multilingual 2.42 74.85 68.90 80.80
103 multilingual 0.16 62.94 59.86 66.02
104 monolingual -1.60 47.42 42.74 52.10
105 multilingual 0.69 65.75 62.33 69.17
106 multilingual 2.52 75.33 69.22 81.45
107 monolingual 0.93 53.27 49.68 56.87
108 monolingual -0.32 50.38 47.18 53.58
109 monolingual -0.80 49.28 45.69 52.86
110 multilingual 1.76 71.37 66.54 76.20
111 multilingual -1.70 53.16 48.63 57.70
112 multilingual -1.45 54.45 50.28 58.63
113 multilingual 2.07 72.99 67.65 78.32
114 monolingual -0.87 49.10 45.42 52.77
115 monolingual -0.60 49.74 46.34 53.13
116 multilingual 1.51 70.02 65.58 74.45
117 multilingual -1.45 54.47 50.30 58.64
118 multilingual -2.33 49.85 44.28 55.41
119 monolingual -0.38 50.24 47.01 53.47
120 multilingual 0.00 62.10 59.05 65.15
121 monolingual -0.31 50.42 47.23 53.61
122 multilingual -0.70 58.40 55.09 61.72
123 monolingual 1.92 55.56 50.57 60.56
124 monolingual -0.14 50.80 47.68 53.92
125 monolingual 0.25 51.71 48.59 54.83
126 multilingual 2.40 74.75 68.84 80.67
127 multilingual 1.46 69.77 65.41 74.13
128 monolingual -1.37 47.95 43.62 52.28
129 monolingual -0.87 49.11 45.44 52.78
130 monolingual 2.44 56.78 50.87 62.69
131 multilingual 0.51 64.78 61.51 68.04
132 monolingual 1.39 54.34 50.16 58.52
133 multilingual -0.55 59.19 55.98 62.40
134 multilingual 1.43 69.63 65.31 73.95
135 multilingual 0.27 63.53 60.41 66.66
136 multilingual 2.16 73.44 67.96 78.92
137 multilingual -1.49 54.26 50.03 58.49
138 multilingual -1.00 56.81 53.20 60.42
139 multilingual -1.94 51.88 46.96 56.80
140 monolingual -1.36 47.97 43.65 52.29
141 multilingual 2.47 75.10 69.07 81.13
142 multilingual -0.93 57.18 53.64 60.71
143 multilingual 1.22 68.49 64.47 72.51
144 monolingual 1.51 54.63 50.27 58.99
145 multilingual -1.89 52.15 47.32 56.99
146 multilingual -1.22 55.69 51.82 59.55
147 multilingual -1.06 56.50 52.82 60.18
148 monolingual -1.91 46.70 41.51 51.89
149 multilingual -1.83 52.48 47.74 57.22
150 multilingual 2.54 75.47 69.32 81.63
151 monolingual 2.52 56.95 50.91 62.99
152 multilingual -1.73 52.99 48.40 57.58
153 monolingual 2.11 56.01 50.69 61.34
154 multilingual 0.47 64.55 61.32 67.78
155 monolingual -0.44 50.11 46.84 53.38
156 monolingual -0.17 50.74 47.61 53.87
157 monolingual 1.12 53.71 49.90 57.53
158 monolingual -2.00 46.49 41.13 51.84
159 monolingual -0.72 49.46 45.95 52.97
160 monolingual 0.99 53.41 49.75 57.08
161 monolingual -0.83 49.20 45.58 52.83
162 multilingual 1.74 71.27 66.47 76.07
163 multilingual -1.34 55.04 51.02 59.07
164 multilingual 0.08 62.49 59.43 65.55
165 monolingual -1.03 48.73 44.86 52.59
166 multilingual -1.45 54.43 50.25 58.62
167 monolingual 2.34 56.54 50.81 62.26
168 monolingual -0.81 49.26 45.66 52.85
169 multilingual -0.02 61.97 58.93 65.02
170 monolingual -2.27 45.86 40.02 51.69
171 multilingual 0.32 63.79 60.65 66.94
172 multilingual 0.81 66.34 62.81 69.87
173 multilingual 0.57 65.08 61.77 68.38
174 multilingual -0.80 57.85 54.44 61.26
175 monolingual 2.04 55.85 50.64 61.05
176 monolingual 0.72 52.78 49.39 56.17
177 monolingual -0.90 49.04 45.33 52.74
178 multilingual -0.47 59.60 56.44 62.76
179 multilingual -1.61 53.60 49.19 58.02
180 multilingual 1.90 72.09 67.04 77.14
181 multilingual 2.35 74.47 68.65 80.29
182 monolingual 0.40 52.06 48.88 55.24
183 monolingual -0.77 49.35 45.79 52.90
184 monolingual 2.57 57.07 50.93 63.21
185 monolingual -1.24 48.25 44.11 52.39
186 multilingual 0.65 65.51 62.13 68.88
187 multilingual -1.87 52.23 47.41 57.04
188 monolingual 0.02 51.17 48.07 54.27
189 monolingual -1.92 46.68 41.48 51.89
190 multilingual -1.61 53.62 49.21 58.03
191 multilingual -1.00 56.83 53.22 60.44
192 monolingual 0.50 52.29 49.06 55.52
193 monolingual 1.24 54.00 50.03 57.98
194 multilingual -1.59 53.74 49.36 58.11
195 multilingual 1.92 72.19 67.11 77.27
196 monolingual 1.13 53.74 49.91 57.56
197 monolingual 1.39 54.34 50.16 58.51
198 multilingual -1.68 53.25 48.74 57.76
199 monolingual -0.62 49.68 46.26 53.10
200 multilingual 0.95 67.10 63.42 70.79
201 multilingual 0.20 63.17 60.08 66.26
202 multilingual -0.67 58.59 55.30 61.88
203 multilingual -1.21 55.71 51.85 59.57
204 monolingual -2.12 46.21 40.64 51.77
205 multilingual -1.23 55.61 51.72 59.49
206 multilingual 2.04 72.81 67.53 78.09
207 multilingual 1.64 70.75 66.11 75.40
208 monolingual 1.32 54.19 50.10 58.27
209 multilingual -1.64 53.46 49.00 57.91
210 multilingual -1.79 52.66 47.98 57.35
211 multilingual 2.46 75.04 69.03 81.05
212 monolingual -0.23 50.58 47.42 53.74
213 monolingual -0.09 50.90 47.79 54.02
214 multilingual -1.59 53.73 49.35 58.11
215 multilingual 0.51 64.78 61.51 68.04
216 multilingual -1.06 56.51 52.83 60.18
217 multilingual -1.26 55.44 51.51 59.36
218 monolingual 1.04 53.53 49.81 57.26
219 multilingual -1.00 56.82 53.22 60.43
220 monolingual 1.64 54.91 50.37 59.46
221 multilingual -1.95 51.85 46.92 56.78
222 monolingual 1.70 55.05 50.41 59.69
223 monolingual -0.28 50.48 47.30 53.66
224 monolingual 1.36 54.28 50.14 58.43
225 multilingual 0.90 66.82 63.19 70.44
226 multilingual -0.19 61.08 58.02 64.14
227 monolingual 0.72 52.79 49.40 56.19
228 monolingual -2.41 45.54 39.46 51.62
229 multilingual -1.33 55.10 51.09 59.11
230 monolingual 1.11 53.69 49.89 57.50
231 monolingual -1.34 48.02 43.74 52.31
232 multilingual 1.66 70.81 66.14 75.47
233 multilingual -0.88 57.48 54.00 60.96
234 monolingual 1.02 53.49 49.79 57.20
235 multilingual 2.25 73.93 68.29 79.58
236 multilingual -1.84 52.43 47.67 57.18
237 multilingual -1.78 52.74 48.08 57.40
238 multilingual 0.98 67.23 63.52 70.95
239 multilingual -0.27 60.67 57.60 63.75
240 monolingual 1.76 55.19 50.46 59.93
241 monolingual 2.44 56.78 50.87 62.69
242 multilingual -2.06 51.23 46.11 56.36
243 monolingual -0.12 50.85 47.74 53.97
244 multilingual 1.09 67.85 63.99 71.71
245 monolingual -1.98 46.54 41.22 51.85
246 monolingual 2.55 57.03 50.92 63.13
247 monolingual -1.15 48.46 44.44 52.48
248 multilingual -2.17 50.68 45.39 55.98
249 multilingual 1.02 67.45 63.68 71.21
250 multilingual 1.52 70.10 65.64 74.55

From predictions to comparisons

comparisons <- comparisons(uncentered_model, variables = "multilingual")

comparisons |> 
    as.data.frame() |> 
    dplyr::select(rowid, contrast, age, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)
rowid contrast age estimate conf.low conf.high
1 multilingual - monolingual 8.29 7.64 2.16 13.12
2 multilingual - monolingual 8.11 7.12 1.31 12.93
3 multilingual - monolingual 9.97 12.59 7.96 17.23
4 multilingual - monolingual 8.34 7.79 2.40 13.18
5 multilingual - monolingual 9.66 11.68 7.28 16.08
6 multilingual - monolingual 10.93 15.43 9.21 21.65
7 multilingual - monolingual 7.84 6.33 -0.04 12.69
8 multilingual - monolingual 9.02 9.81 5.31 14.31
9 multilingual - monolingual 9.36 10.80 6.45 15.15
10 multilingual - monolingual 11.34 16.65 9.50 23.80
11 multilingual - monolingual 11.63 17.50 9.66 25.34
12 multilingual - monolingual 11.41 16.85 9.54 24.17
13 multilingual - monolingual 10.37 13.79 8.61 18.97
14 multilingual - monolingual 11.75 17.86 9.72 26.00
15 multilingual - monolingual 9.58 11.46 7.09 15.84
16 multilingual - monolingual 9.88 12.35 7.79 16.91
17 multilingual - monolingual 8.68 8.81 3.95 13.66
18 multilingual - monolingual 8.74 8.97 4.18 13.76
19 multilingual - monolingual 7.10 4.14 -3.96 12.24
20 multilingual - monolingual 9.51 11.26 6.91 15.62
21 multilingual - monolingual 11.36 16.69 9.51 23.88
22 multilingual - monolingual 7.03 3.94 -4.33 12.21
23 multilingual - monolingual 7.36 4.91 -2.56 12.37
24 multilingual - monolingual 7.82 6.27 -0.14 12.68
25 multilingual - monolingual 10.85 15.21 9.14 21.27
26 multilingual - monolingual 10.68 14.69 8.97 20.41
27 multilingual - monolingual 11.86 18.18 9.76 26.59
28 multilingual - monolingual 9.33 10.73 6.37 15.08
29 multilingual - monolingual 7.37 4.94 -2.49 12.38
30 multilingual - monolingual 10.24 13.42 8.43 18.40
31 multilingual - monolingual 10.79 15.03 9.09 20.98
32 multilingual - monolingual 7.69 5.87 -0.84 12.58
33 multilingual - monolingual 8.98 9.70 5.17 14.23
34 multilingual - monolingual 8.12 7.17 1.38 12.95
35 multilingual - monolingual 7.29 4.70 -2.93 12.34
36 multilingual - monolingual 8.98 9.69 5.15 14.22
37 multilingual - monolingual 7.32 4.80 -2.75 12.36
38 multilingual - monolingual 8.13 7.18 1.41 12.95
39 multilingual - monolingual 7.27 4.65 -3.02 12.33
40 multilingual - monolingual 10.35 13.73 8.59 18.88
41 multilingual - monolingual 8.49 8.24 3.10 13.38
42 multilingual - monolingual 7.50 5.33 -1.79 12.46
43 multilingual - monolingual 7.36 4.91 -2.56 12.37
44 multilingual - monolingual 11.40 16.83 9.54 24.12
45 multilingual - monolingual 10.77 14.97 9.07 20.87
46 multilingual - monolingual 11.08 15.89 9.33 22.45
47 multilingual - monolingual 11.91 18.33 9.79 26.88
48 multilingual - monolingual 7.52 5.37 -1.72 12.47
49 multilingual - monolingual 7.50 5.31 -1.84 12.45
50 multilingual - monolingual 10.99 15.63 9.26 21.99
51 multilingual - monolingual 10.92 15.42 9.20 21.63
52 multilingual - monolingual 7.05 3.99 -4.25 12.22
53 multilingual - monolingual 10.90 15.34 9.18 21.49
54 multilingual - monolingual 10.65 14.60 8.94 20.26
55 multilingual - monolingual 10.15 13.14 8.29 17.99
56 multilingual - monolingual 9.40 10.94 6.59 15.28
57 multilingual - monolingual 7.78 6.16 -0.33 12.65
58 multilingual - monolingual 7.04 3.97 -4.28 12.22
59 multilingual - monolingual 9.26 10.52 6.15 14.89
60 multilingual - monolingual 9.46 11.11 6.76 15.45
61 multilingual - monolingual 8.95 9.59 5.03 14.15
62 multilingual - monolingual 9.32 10.70 6.35 15.05
63 multilingual - monolingual 10.57 14.37 8.85 19.88
64 multilingual - monolingual 7.28 4.66 -3.00 12.33
65 multilingual - monolingual 8.77 9.08 4.34 13.82
66 multilingual - monolingual 11.01 15.69 9.28 22.10
67 multilingual - monolingual 11.18 16.17 9.40 22.95
68 multilingual - monolingual 8.19 7.35 1.69 13.01
69 multilingual - monolingual 8.77 9.07 4.32 13.81
70 multilingual - monolingual 11.28 16.48 9.47 23.50
71 multilingual - monolingual 11.27 16.44 9.46 23.42
72 multilingual - monolingual 8.48 8.21 3.06 13.36
73 multilingual - monolingual 7.74 6.02 -0.58 12.61
74 multilingual - monolingual 10.52 14.23 8.80 19.66
75 multilingual - monolingual 7.52 5.38 -1.71 12.47
76 multilingual - monolingual 7.17 4.34 -3.59 12.28
77 multilingual - monolingual 12.00 18.59 9.82 27.35
78 multilingual - monolingual 7.17 4.36 -3.56 12.28
79 multilingual - monolingual 8.69 8.84 3.99 13.68
80 multilingual - monolingual 11.58 17.34 9.63 25.05
81 multilingual - monolingual 10.09 12.95 8.18 17.72
82 multilingual - monolingual 8.43 8.07 2.84 13.30
83 multilingual - monolingual 10.69 14.73 8.99 20.47
84 multilingual - monolingual 11.17 16.15 9.39 22.90
85 multilingual - monolingual 8.57 8.48 3.47 13.49
86 multilingual - monolingual 9.46 11.11 6.76 15.46
87 multilingual - monolingual 10.49 14.13 8.76 19.50
88 multilingual - monolingual 10.21 13.31 8.38 18.24
89 multilingual - monolingual 10.22 13.34 8.40 18.29
90 multilingual - monolingual 11.89 18.27 9.78 26.76
91 multilingual - monolingual 9.07 9.96 5.50 14.43
92 multilingual - monolingual 7.60 5.61 -1.30 12.52
93 multilingual - monolingual 9.63 11.60 7.21 16.00
94 multilingual - monolingual 8.13 7.17 1.39 12.95
95 multilingual - monolingual 9.43 11.02 6.68 15.37
96 multilingual - monolingual 8.85 9.31 4.65 13.96
97 multilingual - monolingual 11.92 18.35 9.79 26.91
98 multilingual - monolingual 8.94 9.57 5.01 14.14
99 multilingual - monolingual 8.15 7.23 1.49 12.97
100 multilingual - monolingual 10.12 13.04 8.23 17.85
101 multilingual - monolingual 7.68 5.86 -0.85 12.58
102 multilingual - monolingual 11.84 18.11 9.76 26.47
103 multilingual - monolingual 9.58 11.44 7.07 15.81
104 multilingual - monolingual 7.82 6.25 -0.17 12.67
105 multilingual - monolingual 10.11 13.02 8.22 17.82
106 multilingual - monolingual 11.93 18.39 9.79 26.98
107 multilingual - monolingual 10.34 13.71 8.58 18.84
108 multilingual - monolingual 9.09 10.03 5.58 14.47
109 multilingual - monolingual 8.62 8.62 3.67 13.56
110 multilingual - monolingual 11.18 16.16 9.40 22.93
111 multilingual - monolingual 7.72 5.97 -0.66 12.60
112 multilingual - monolingual 7.96 6.69 0.59 12.80
113 multilingual - monolingual 11.48 17.07 9.58 24.56
114 multilingual - monolingual 8.54 8.39 3.33 13.45
115 multilingual - monolingual 8.82 9.20 4.51 13.90
116 multilingual - monolingual 10.92 15.41 9.20 21.62
117 multilingual - monolingual 7.97 6.70 0.60 12.80
118 multilingual - monolingual 7.09 4.11 -4.02 12.24
119 multilingual - monolingual 9.03 9.84 5.35 14.33
120 multilingual - monolingual 9.42 10.97 6.63 15.32
121 multilingual - monolingual 9.11 10.07 5.63 14.51
122 multilingual - monolingual 8.71 8.90 4.09 13.72
123 multilingual - monolingual 11.33 16.63 9.50 23.75
124 multilingual - monolingual 9.28 10.56 6.19 14.92
125 multilingual - monolingual 9.67 11.72 7.31 16.13
126 multilingual - monolingual 11.82 18.06 9.75 26.38
127 multilingual - monolingual 10.87 15.27 9.16 21.38
128 multilingual - monolingual 8.04 6.93 0.99 12.87
129 multilingual - monolingual 8.54 8.40 3.35 13.45
130 multilingual - monolingual 11.86 18.17 9.76 26.58
131 multilingual - monolingual 9.92 12.47 7.88 17.07
132 multilingual - monolingual 10.80 15.07 9.10 21.04
133 multilingual - monolingual 8.86 9.34 4.70 13.99
134 multilingual - monolingual 10.85 15.19 9.14 21.25
135 multilingual - monolingual 9.69 11.78 7.36 16.19
136 multilingual - monolingual 11.57 17.33 9.63 25.02
137 multilingual - monolingual 7.93 6.58 0.40 12.76
138 multilingual - monolingual 8.41 8.01 2.75 13.27
139 multilingual - monolingual 7.47 5.25 -1.95 12.44
140 multilingual - monolingual 8.05 6.95 1.03 12.88
141 multilingual - monolingual 11.89 18.26 9.78 26.74
142 multilingual - monolingual 8.48 8.22 3.07 13.37
143 multilingual - monolingual 10.63 14.55 8.92 20.18
144 multilingual - monolingual 10.93 15.43 9.21 21.66
145 multilingual - monolingual 7.53 5.40 -1.67 12.47
146 multilingual - monolingual 8.20 7.38 1.74 13.02
147 multilingual - monolingual 8.35 7.84 2.48 13.20
148 multilingual - monolingual 7.51 5.34 -1.78 12.46
149 multilingual - monolingual 7.59 5.59 -1.34 12.51
150 multilingual - monolingual 11.96 18.47 9.80 27.13
151 multilingual - monolingual 11.93 18.39 9.79 26.98
152 multilingual - monolingual 7.69 5.87 -0.84 12.58
153 multilingual - monolingual 11.53 17.20 9.61 24.79
154 multilingual - monolingual 9.88 12.35 7.79 16.90
155 multilingual - monolingual 8.98 9.68 5.14 14.21
156 multilingual - monolingual 9.25 10.48 6.11 14.85
157 multilingual - monolingual 10.53 14.27 8.82 19.72
158 multilingual - monolingual 7.41 5.06 -2.28 12.40
159 multilingual - monolingual 8.70 8.85 4.01 13.69
160 multilingual - monolingual 10.40 13.89 8.66 19.12
161 multilingual - monolingual 8.58 8.52 3.53 13.51
162 multilingual - monolingual 11.16 16.11 9.38 22.84
163 multilingual - monolingual 8.08 7.02 1.14 12.90
164 multilingual - monolingual 9.49 11.19 6.84 15.54
165 multilingual - monolingual 8.38 7.92 2.60 13.23
166 multilingual - monolingual 7.96 6.68 0.57 12.79
167 multilingual - monolingual 11.75 17.87 9.72 26.01
168 multilingual - monolingual 8.61 8.59 3.63 13.55
169 multilingual - monolingual 9.39 10.90 6.56 15.25
170 multilingual - monolingual 7.14 4.26 -3.74 12.26
171 multilingual - monolingual 9.74 11.92 7.47 16.37
172 multilingual - monolingual 10.22 13.35 8.40 18.30
173 multilingual - monolingual 9.98 12.64 7.99 17.29
174 multilingual - monolingual 8.61 8.59 3.64 13.55
175 multilingual - monolingual 11.46 16.99 9.57 24.41
176 multilingual - monolingual 10.13 13.08 8.25 17.91
177 multilingual - monolingual 8.51 8.31 3.22 13.41
178 multilingual - monolingual 8.94 9.57 5.01 14.14
179 multilingual - monolingual 7.80 6.21 -0.24 12.66
180 multilingual - monolingual 11.31 16.57 9.48 23.65
181 multilingual - monolingual 11.77 17.90 9.72 26.08
182 multilingual - monolingual 9.82 12.16 7.66 16.66
183 multilingual - monolingual 8.65 8.71 3.80 13.61
184 multilingual - monolingual 11.98 18.54 9.81 27.27
185 multilingual - monolingual 8.17 7.31 1.62 13.00
186 multilingual - monolingual 10.06 12.88 8.14 17.63
187 multilingual - monolingual 7.54 5.44 -1.60 12.48
188 multilingual - monolingual 9.44 11.03 6.68 15.38
189 multilingual - monolingual 7.50 5.31 -1.83 12.45
190 multilingual - monolingual 7.81 6.22 -0.22 12.67
191 multilingual - monolingual 8.41 8.02 2.76 13.28
192 multilingual - monolingual 9.92 12.46 7.87 17.05
193 multilingual - monolingual 10.66 14.64 8.95 20.32
194 multilingual - monolingual 7.83 6.29 -0.11 12.68
195 multilingual - monolingual 11.33 16.63 9.50 23.75
196 multilingual - monolingual 10.54 14.30 8.83 19.77
197 multilingual - monolingual 10.80 15.06 9.10 21.03
198 multilingual - monolingual 7.74 6.02 -0.58 12.61
199 multilingual - monolingual 8.79 9.13 4.40 13.85
200 multilingual - monolingual 10.37 13.78 8.61 18.95
201 multilingual - monolingual 9.62 11.57 7.19 15.96
202 multilingual - monolingual 8.75 9.01 4.23 13.78
203 multilingual - monolingual 8.20 7.39 1.76 13.03
204 multilingual - monolingual 7.29 4.71 -2.93 12.34
205 multilingual - monolingual 8.18 7.34 1.67 13.01
206 multilingual - monolingual 11.45 16.97 9.57 24.38
207 multilingual - monolingual 11.06 15.82 9.31 22.33
208 multilingual - monolingual 10.74 14.87 9.03 20.71
209 multilingual - monolingual 7.77 6.13 -0.38 12.64
210 multilingual - monolingual 7.62 5.69 -1.16 12.53
211 multilingual - monolingual 11.87 18.22 9.77 26.67
212 multilingual - monolingual 9.18 10.28 5.88 14.68
213 multilingual - monolingual 9.32 10.69 6.34 15.04
214 multilingual - monolingual 7.83 6.28 -0.11 12.68
215 multilingual - monolingual 9.92 12.47 7.88 17.07
216 multilingual - monolingual 8.35 7.84 2.48 13.20
217 multilingual - monolingual 8.15 7.24 1.51 12.97
218 multilingual - monolingual 10.46 14.04 8.72 19.36
219 multilingual - monolingual 8.41 8.02 2.76 13.28
220 multilingual - monolingual 11.05 15.80 9.31 22.29
221 multilingual - monolingual 7.47 5.23 -1.97 12.44
222 multilingual - monolingual 11.11 15.97 9.35 22.59
223 multilingual - monolingual 9.14 10.15 5.73 14.57
224 multilingual - monolingual 10.78 14.99 9.08 20.91
225 multilingual - monolingual 10.31 13.62 8.53 18.70
226 multilingual - monolingual 9.22 10.40 6.02 14.79
227 multilingual - monolingual 10.14 13.10 8.26 17.93
228 multilingual - monolingual 7.00 3.85 -4.49 12.20
229 multilingual - monolingual 8.09 7.05 1.19 12.91
230 multilingual - monolingual 10.52 14.24 8.81 19.68
231 multilingual - monolingual 8.08 7.02 1.14 12.90
232 multilingual - monolingual 11.07 15.85 9.32 22.38
233 multilingual - monolingual 8.54 8.39 3.33 13.45
234 multilingual - monolingual 10.44 13.99 8.70 19.28
235 multilingual - monolingual 11.66 17.60 9.68 25.53
236 multilingual - monolingual 7.58 5.55 -1.40 12.51
237 multilingual - monolingual 7.64 5.73 -1.08 12.54
238 multilingual - monolingual 10.39 13.85 8.64 19.06
239 multilingual - monolingual 9.14 10.17 5.75 14.59
240 multilingual - monolingual 11.17 16.15 9.39 22.91
241 multilingual - monolingual 11.86 18.17 9.76 26.58
242 multilingual - monolingual 7.35 4.89 -2.60 12.37
243 multilingual - monolingual 9.30 10.63 6.27 14.99
244 multilingual - monolingual 10.51 14.19 8.79 19.60
245 multilingual - monolingual 7.43 5.13 -2.16 12.42
246 multilingual - monolingual 11.96 18.49 9.81 27.17
247 multilingual - monolingual 8.27 7.58 2.06 13.10
248 multilingual - monolingual 7.25 4.58 -3.16 12.32
249 multilingual - monolingual 10.43 13.97 8.69 19.25
250 multilingual - monolingual 10.93 15.45 9.21 21.69

Calcualting average marginal effects

AME <- avg_comparisons(uncentered_model)

AME |> 
    as.data.frame() |> 
    dplyr::select(term, contrast, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)
term contrast estimate conf.low conf.high
age mean(+1) 3.81 2.34 5.29
multilingual mean(multilingual) - mean(monolingual) 10.97 6.62 15.31

marginaleffects can also make plots

plot_predictions(uncentered_model, condition = list("multilingual"))

marginaleffects can also make plots

plot_predictions(uncentered_model, condition = list("age"))

marginaleffects can also make plots

plot_predictions(uncentered_model, condition = list("age", "multilingual"))

marginaleffects can also make plots

Because the plots are made with ggplot2, we can customize them as we like.

plot_predictions(uncentered_model, condition = list("age", "multilingual")) +
    theme(legend.title = element_blank()) +
    scale_color_manual(values = c("monolingual" = "royalblue", 
                                "multilingual" = "darkorange")) +
    scale_fill_manual(values = c("monolingual" = "royalblue", 
                                "multilingual" = "darkorange")) +
    labs(x = "Age", y = "Cognitive flexibility score")

Generalized linear models

  • GLMs such as logistic regression and Poisson regression are also supported by marginaleffects.

  • Good news for us, because everything interacts in a GLM due to the link function (e.g., logit, log, softmax).

    • Extremely easy to mislead ourselves if we try to interpret parameters.

Logistic regression

  • Hypothetical dataset: Touchscreen experiment with Chimpanzees
    • After viewing a stimulus, the chimpanzee touches the screen on either the ‘Agent’ or the ‘Patient’ in a scence.
    • Estimand: what is the average difference in the probability of touching the ‘Agent’ side between two conditions (‘experimental’ vs ‘control’)?
      • \(\frac{1}{N} \sum_{i=1}^N [P(Y_i = \text{Agent} | \text{experimental}) - P(Y_i = \text{Agent} | \text{control})]\)

Logistic regression

\[ \text{Agent} \sim \text{Bernoulli}(\pi) \] \[ \pi = \text{logit}^{-1}(\alpha + \beta \times \text{experimental}) \]


Call:
glm(formula = agent ~ condition, family = binomial(link = "logit"), 
    data = d)

Coefficients:
                      Estimate Std. Error z value Pr(>|z|)  
(Intercept)             0.5465     0.3789   1.443   0.1491  
conditionexperimental   1.3253     0.6573   2.016   0.0438 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 67.48  on 59  degrees of freedom
Residual deviance: 62.99  on 58  degrees of freedom
AIC: 66.99

Number of Fisher Scoring iterations: 4

Logistic regression

  • Even with just a single predictor, we cannot get our estimate from looking at model coefficients.

  • Due to the link function, the effect of the experimental condition depends on the intercept.

  • ‘Mechanistic’ rather than conceptual interaction.

    • To get the AME we need to calculate:
      • \(\frac{1}{1 + e^{-(\alpha + \beta \times \text{experimental})}} - \frac{1}{1 + e^{-\alpha}}\)

Logistic regression

The marginal effect on the probability scale for a log odds ratio of 1.33 depends strongly on the intercept!

marginaleffects to the rescue

AME <- avg_comparisons(logistic_model, variables = "condition")

AME |> 
    as.data.frame() |> 
    dplyr::select(term, contrast, estimate, conf.low, conf.high) |> 
    print(digits = 2)
       term                           contrast estimate conf.low conf.high
1 condition mean(experimental) - mean(control)     0.23    0.022      0.44
plot_predictions(logistic_model, condition = list("condition")) + ylab("Pr(Agent)")

  • NOTE: I should write out the estimate in line with the original goal to make the outputs easier to report and interpret.

But how does it work?

  • For frequentist models, marginaleffects by default uses the delta method (Taylor series expansion)

    • Various bootstrap and weighting methods supported
  • For Bayesian models (e.g., fit with brms), marginaleffects, uses full posterior draws

  • More info: https://marginaleffects.com/chapters/uncertainty.html

Review thus far

Reporting model results using marginaleffects:

  • Makes your model results easier to communicate
  • Avoid misinterpreting the magnitude of effects
  • Avoid misinterpreting the direction of effects

Non-linear effects

Example: relationship between population density and structural diversity in language (measured by normalized entropy)

Non-linear effects


Family: gaussian 
Link function: identity 

Formula:
structural_diversity ~ s(log_density, k = 8)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0.676467   0.009797   69.05   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Approximate significance of smooth terms:
                 edf Ref.df     F p-value    
s(log_density) 6.221  6.798 49.08  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

R-sq.(adj) =  0.623   Deviance explained = 63.5%
GCV = 0.019914  Scale est. = 0.019195  n = 200

Non-linear effects

slopes_df <- slopes(gam_model) |> 
    as.data.frame() |> 
    dplyr::select(estimate, conf.low, conf.high, log_density, structural_diversity)

print(head(slopes_df), digits = 2)
  estimate conf.low conf.high log_density structural_diversity
1   0.1590    0.110     0.208        0.88                 0.45
2   0.0023   -0.059     0.064        5.88                 0.89
3  -0.0065   -0.059     0.046        2.09                 0.68
4   0.0392   -0.020     0.098        6.83                 0.83
5   0.0912   -0.015     0.197        7.40                 0.96
6   0.0240   -0.100     0.148       -1.54                 0.37

Non-linear effects: average slopes

\(\mathbb{E}\left[\frac{\partial \text{structural diversity}}{\partial \text{ log pop density}}\right]\)

avg_slopes(gam_model) |> 
    as.data.frame() |> 
    dplyr::select(term, estimate, conf.low, conf.high) |> 
    knitr::kable(digits = 2)
term estimate conf.low conf.high
log_density 0.07 0.06 0.08

The bigger picture

  • Directly interpreting statistical models is a full of traps that can mislead us. Working with predictions and marginal effects allows us to avoid many of these traps.

  • marginaleffects can be used with a huge array of models:

    • linear models, generalized linear models (e.g., lm, glm)
    • generalized additive models (e.g., mgcv, gamm4)
    • multilevel models (e.g., lme4, glmmTMB)
    • Bayesian models (e.g., brms)
    • machine learning models (e.g., tidymodels, mlr3)

Can marginal effects improve the reliability of science?

  • On the one hand, we can avoid falling into many of the traps of interpreting model parameters.

  • On the other hand, this way of doing things brings the estimand into focus. We are being more explicit about the quantity we are trying to estimate–which is not tied to any particular model.

Red Card Dataset

2,053 soccer players in the first male divisions of England, Germany, France, and Spain in the 2012–2013 season.

Outcome: number of red cards recieved for 146,028 player-referee dyads.

Covariates: player and referee traits, match statistics

Many analysts, many estimates

“Are soccer players with dark skin tone more likely than those with light skin tone to receive red cards from referees?”

Many analysts, many estimates

“Are soccer players with dark skin tone more likely than those with light skin tone to receive red cards from referees?”

Who watches the watchmen? Reanalysis of the meta-analysis

Many analysts, many estimands

“Thus, we argue that the results obtained in the CSI showed such a large variation because the 29 teams pursued (at least) four different research questions and therefore used different research designs.”

Controlled direct effect: \(E[Y(\text{dark}) - Y(\text{light}) | M = m]\); 486 parametric models

Causal (not statistical) assumptions needed to inform research design

Three critical choices in quantitative research. From: Lundberg, I., Johnson, R., & Stewart, B. M. (2021). What is your estimand? Defining the target quantity connects statistical evidence to theory. American Sociological Review, 86(3), 532-565.

Breakout session (10-15 minutes)

Goals:

  1. Write down the estimand for one of your research questions
    • Remember: this is a unit-specific comparison, ideally with reference to a specific population.
  2. Draft a statistical model that you could use as your estimator
    • Write down the marginaleffects code that would yield your estimate
    • If you don’t have the relevant software on your laptop, find a neighbor who does.
    • Pseudocode is fine!

Report-back/Q&A session

  • How did it go? What friction did you encounter?

Part 2: Towards a less stressful computational workflow

Statistics can feel like a thankless list of things that all have to be right or else nothing is right. - Not totally wrong, but leads to anxiety

  • Calls to improve statistical practice in science imply time commitments:
    • More time for training
    • More time for coding tasks (recall the multiverse)

Academic time allocation

Brauer, K. (2021). ” I’ll Finish It This Week” And Other Lies. arXiv preprint arXiv:2103.16574.

Harsh realities

  • Time allocation is a zero-sum game
    • Time spent on doing better statistics is time that could be spent on the many other aspects of science and professional development
  • If raise the bar on statistical pracitce we need workflows that are reproducible, understandable, and which reduce cognitive burden.

The tangled web of data analysis

  • What happens when Reviewer #2 asks for a change to the data/model?

targets

Package for reproducible research pipelines by Will Landau:

  • Encodes every step of your data analysis as a directed acyclic graph

  • Automatically detection of dependencies, which means:

    • If you change some part of your pipeline, targets will rerun (only) the parts that depend on the changed part.

Cognitive flexibility example

# _targets.R file
tar_source(files = "R") # load the R scripts in my dir
tar_option_set(packages = c("tidyverse", "marginaleffects")) # load any packages used
# define pipeline
list(
    tar_target(data_file, "data/cognitive_flexibility.csv", format = "file"),
    tar_target(data, read_csv(data_file)),
    tar_target(model, fit_model(data)),
    tar_target(marginal_effect, avg_comparisons(model, newdata = data, variables = "multilingual")),
    tar_target(model_plot, plot_predictions(model, condition = list("multilingual")))
)

tar_visnetwork()

tar_visnetwork()

Looking inside in the pipeline?

# fit_model.R
fit_model <- function(data, age_moderator = TRUE){
    if(age_moderator){
        model <- lm(cognitive_flexibility ~ multilingual * age, data = data)
    } else {
        model <- lm(cognitive_flexibility ~ multilingual, data = data)
    }
    return(model)
}

Running the pipeline

tar_make()
▶ dispatched target data_file
● completed target data_file [0.402 seconds, 17.112 kilobytes]
▶ dispatched target data
Rows: 250 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): multilingual
dbl (3): cognitive_flexibility, age, age_c

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
● completed target data [0.114 seconds, 5.687 kilobytes]
▶ dispatched target model
● completed target model [0.003 seconds, 17.378 kilobytes]
▶ dispatched target marginal_effect
● completed target marginal_effect [0.089 seconds, 41.914 kilobytes]
▶ dispatched target model_plot
● completed target model_plot [0.036 seconds, 253.351 kilobytes]
▶ ended pipeline [0.809 seconds]

Re-inspect the pipeline

tar_visnetwork()

See what happens if we change the data

# remove anomalous observations
problematic_obs <- c(3, 15, 27)

write_csv(data[-problematic_obs,], "data/cognitive_flexibility.csv")
tar_visnetwork()

See what happens if we change the data

tar_outdated()
[1] "marginal_effect" "data"            "model_plot"      "data_file"      
[5] "model"          

See what happens if we change the model function

# R/fit_model.R, now with quadratic effect for age
fit_model <- function(data, age_moderator = TRUE){
    if(age_moderator){
        model <- lm(cognitive_flexibility ~ multilingual * poly(age, 2), data = data)
    } else {
        model <- lm(cognitive_flexibility ~ multilingual, data = data)
    }
    return(model)
}

See what happens if we change the model function

tar_visnetwork()

See what happens if we change the model function

tar_outdated()
[1] "marginal_effect" "model_plot"      "model"          

Real world example

  • Big latent, multilevel timeseries model (>10,000 parameters)
  • Many outputs (effect sizes, constrasts, model visualizations, data dashboard)
  • Entire pipeline takes ~6 days to run

Targets at scale

Targets at scale

library(targets)
source("R/functions.R")

# Dependencies for this pipeline
tar_option_set(packages = c(
  "tidyverse",
  "rethinking",
  "rstan",
  "patchwork",
  "ggsci",
  "RColorBrewer",
  "cubelyr",
  "proxy",
  "foreach",
  "doParallel",
  "ggdist"
))

list(
  # Extract raw data
  tar_target(rdata_file,
             "data/input-to-models.RData", format = "file"),
  tar_target(footage_file,
             "data/Study1_footage_details.xlsx"),
  tar_target(infants_file,
             "data/Infant_age.xlsx"),
  tar_target(d_raw_time,
             extract_raw_data(rdata_file, footage_file, infants_file, "time")),
  
  ## Generate data for Shiny dashboard #####
  tar_target(d_dash,
             dashboard_data(d_raw_time)),

  # Data for time series model 
  tar_target(d_time,
             time_data(d_raw_time)),
  
  ## Generate data for time-averaged model comparison
  tar_target(d_avg,
             time_data_avg(d_time)),
  
  ###### Time Averaged Analyses #########################
  tar_target(model_compare,
             avg_model_comparison(d_avg)),
  
  ###### Time Series Analyses ##########################
  # Fit time series model
  tar_target(fit_time,
             fit_time_series(d_time, n_iter = 1000, n_chains = 8)),
  # Make predictions for time series
  tar_target(pred_time,
             pred_species_pid_time(fit_time, d_time)),
  
  # Make plots + export summaries of posterior predictions to .csv
  # time series plots
  tar_target(p_species_time_agent,
             plot_species_time_AOI(pred_time, "Agent")),
  tar_target(p_species_time_patient,
             plot_species_time_AOI(pred_time, "Patient")),
  tar_target(p_species_time_other,
             plot_species_time_AOI(pred_time, "Other")),
  tar_target(p_species_time_OR,
             plot_species_time_OR(pred_time)),
  tar_target(p_species_time_APO,
             plot_species_time_OR_APO(pred_time)),
  tar_target(p_species_time_AO,
             plot_species_time_OR_AO(pred_time)),
  # AOI size plots
  tar_target(p_species_AOI_size,
             plot_species_AOI_size(fit_time, d_time)),
  # AP movement diff
  tar_target(p_species_AP_md,
             plot_species_AP_md(fit_time, d_time))
  
)

Advanced use of targets

  • Dynamic branching patterns, when there are too many targets to specify manually

  • Great support for parallel processing using future

Conclusion / Q&A